getwd()
setwd("/Users/JONGHYUK/Desktop")
library(foreign)
library(arm)
cws <- read.dta("final1.dta")

######Regime Classification
##1. Social democracy (Solid-lines): 5, 6, 12, 14, 15    
##2. Liberal(Dash-lines):  1, 4, 9, 11, 13, 16, 17, 18
##3. Conservative(Dotted-lines): 2, 3, 7, 8, 10

##Multiple Inputations by Amerlia II (http://gking.harvard.edu/amelia/)

require(Amelia)
second.data <- amelia(second,ts="year",cs="idn")

head(second.data)
save(second.data, file="final.RData")
write.amelia(obj=second.data, file.stem = "final", format="dta")

missmap(second.data)
plot(second.data, which.var=24)
plot(second.data)
overimpute(second.data, var="")
disperse(second.data, dims=1, m=5)
tscsPlot
read.table("D:\Document\2010_Gary King\project")
a <- read.table("final.Rdata")
head(final)
ls()

**zelig with ls.mixed
library(Zelig)
a<-read.dta("final1_recode.dta")
z.out1<-zelig(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+tag(1|ltcabcuml), data=a, model="ls.mixed")
summary(z.out1)



***data setting
setwd("D:/Document/2010_Gary King/project")
library(foreign)
a <-read.dta("final1_recode.dta")
b<-read.dta("final2_recode.dta")
c<-read.dta("final3_recode.dta")
d<-read.dta("final4_recode.dta")
e<-read.dta("final5_recode.dta")

na.a <- !is.na( apply (a [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
a.nomiss <- a [na.a ,]
na.b <- !is.na( apply (b [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
b.nomiss <- b [na.b ,]
na.c <- !is.na( apply (c [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
c.nomiss <- c [na.c ,]
na.d <- !is.na( apply (d [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
d.nomiss <- d [na.d ,]
na.e <- !is.na( apply (e [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
e.nomiss <- e [na.e ,]


**lme

library(nlme)
jing.a<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=a.nomiss)
jing.b<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=b.nomiss)
jing.c<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=c.nomiss)
jing.d<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=d.nomiss)
jing.e<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=e.nomiss)

##other dvs

ue.bothl<-lme(uescorecov ~ ltcabcuml+cncabcuml+wmobcal + vturnl+ oldl+lgdppcl+unempl+dinvocl+openkl, data=jing.a, random=~ltcabcuml+cncabcuml)
ue.botha<-update(ue.bothl, correlation=corAR1())
summary(ue.botha)

ue.bothl<-lme(decom ~ ltcabcuml+cncabcuml, data=jing.a, random=~ltcabcuml+cncabcuml)


na.a <- !is.na( apply (a [,c("decom", "ue", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
a.nomiss <- a [na.a ,]
jing.a<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=a.nomiss)
gross.both<-lme(gross~ltcabcuml+cncabcuml+veto1l+wmobcal+oldl+lgdppcl+unempl+dinvocl+openkl, data=jing.a, random=~ltcabcuml+cncabcuml)
summary(gross.both)
auto.gross<-update(gross.both, correlation=corAR1())
summary(
lme.gross<-lme(gross~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl, data=jing.a, random=~ltcabcuml+cncabcuml)


ue.bothl<-lme(uewait ~ ltcabcuml + cncabcuml + wmobcal + vturnl + oldl + cpil +mill + dinvocl + openkl, data=jing.a, random=~ltcabcuml+cncabcuml)
ue.botha<-update(ue.bothl, correlation=corAR1())
summary(ue.botha)

##gross-test
gross.intf<-lmer(formula = gross ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl +openkl+(1 | idn), data = a.nomiss)
gross.ints<-lmer(formula = gross ~ ltcabcuml + cncabcuml + vturnl + oldl + lgdppcl+ cpil + mill + dinvocl +openkl+(1 | idn), data = a.nomiss)
summary(gross.intf)
boths<-lmer(formula = uescorecov ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl +openkl+(1 + ltcabcuml | idn)+(1 + cncabcuml | idn), data = a.nomiss)
gross.bothl<-lme(gross~ltcabcuml + cncabcuml + wmobcal + vturnl + oldl + lgdppcl + unempl, data=jing.a, random=~ltcabcuml+cncabcuml)
gross.botha<-update(gross.bothl, correlation=corAR1())
summary(boths)

ue.intf<-lmer(formula = ue ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl +openkl+(1 | idn), data = a.nomiss)
ue.ints<-lmer(formula = ue ~ ltcabcuml + cncabcuml + vturnl + oldl + lgdppcl+cpil + mill + dinvocl +openkl+(1 | idn), data = a.nomiss)
ue.boths<-lmer(formula = ue ~ ltcabcuml + cncabcuml + vturnl + oldl + cpil + mill + dinvocl +openkl+(1 + ltcabcuml | idn)+(1 + cncabcuml | idn), data = jing.a)
ue.bothl<-lme(ue ~ ltcabcuml + cncabcuml + vturnl + oldl + cpil + mill + dinvocl + openkl, data=jing.a, random=~ltcabcuml+cncabcuml)
ue.botha<-update(ue.bothl, correlation=corAR1())
summary(ue.intf)

summary(gross.botha)
summary(gross.int)


lme.a<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl, data=jing.a, random=~ltcabcuml+cncabcuml)
lmer.a<-lmer(formula = decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 + ltcabcuml | idn), data = a.nomiss)
summary(lmer.a)
summary(lme.a)


###lmer

library(lme4)
lmer.a<-lmer(formula = decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 + ltcabcuml | idn) + (1 + cncabcuml | idn), data = a.nomiss)
summary(lmer.a)
summary(lme.a)
gross.int<-lmer(formula = gross ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl +openkl+(1 | idn), data = a.nomiss)
summary(gross.int)
gross.both<-lmer(formula = gross ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 + ltcabcuml | idn) + (1 + cncabcuml | idn), data = a.nomiss)
summary(gross.both)
lme.gross<-lme(gross~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl, data=jing.a, random=~ltcabcuml+cncabcuml)

ue.int<-lmer(formula = ue ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 | idn), data = a.nomiss)
summary(ue.int)
ue.both<-lmer(formula = ue ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 + ltcabcuml | idn) + (1 + cncabcuml | idn), data = a.nomiss)
summary(ue.both)
coef(ue.both)
head(a.nomiss)

act.int<-lmer(formula = act1_10 ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 | idn), data = a.nomiss)
summary(act.int)
act.both<-lmer(formula = act1_10 ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 + ltcabcuml | idn) + (1 + cncabcuml | idn), data = a.nomiss)
summary(act.both)

almp.int<-lmer(formula = almp ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 | idn), data = a.nomiss)
summary(almp.int)
almp.both<-lmer(formula = almp ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 + ltcabcuml | idn) + (1 + cncabcuml | idn), data = a.nomiss)
summary(almp.both)

uewait.int<-lmer(formula = uewait ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 | idn), data = a.nomiss)
summary(uewait.int)
uewait.both<-lmer(formula = uewait ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 + ltcabcuml | idn) + (1 + cncabcuml | idn), data = a.nomiss)
summary(uewait.both)

uedurc.int<-lmer(formula = uedurc ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 | idn), data = a.nomiss)
summary(uedurc.int)
uedurc.both<-lmer(formula = uedurc ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 + ltcabcuml | idn) + (1 + cncabcuml | idn), data = a.nomiss)
summary(uedurc.both)

taxpay.int<-lmer(formula = taxpay~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 | idn), data = a.nomiss)
summary(taxpay.int)
taxpay.both<-lmer(formula = taxpay ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 + ltcabcuml | idn) + (1 + cncabcuml | idn), data = a.nomiss)
summary(taxpay.both)

epl.int<-lmer(formula = epl~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 | idn), data = a.nomiss)
summary(epl.int)
epl.both<-lmer(formula = epl ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl +(1 + ltcabcuml | idn) + (1 + cncabcuml | idn), data = a.nomiss)
summary(epl.both)

jing.a
both.a<-(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn+year, data=a.nomiss)
ergoStool
ergoStool.mat<-asTable(ergoStool)
ergoStool.mat
ergoStool.new<-balancedGrouped(effort ~ Type|Subject, data=ergoStool.mat)

table(a.nomiss$idn)
table(b.nomiss$idn)
table(getGroups(b.nomiss))
unique(table(a.nomiss$idn))
table(isBalanced(jing.a))
balance.a<-asTable(jing.a)
balance.a<-asTable(jing.a)
names()
test<-read.dta("test.dta")
table(getGroups(test))
names(test)

lme.a<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl, data=jing.a, random=~ltcabcuml)
lme.b<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl, data=jing.b, random=~ltcabcuml)
lme.c<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl, data=jing.c, random=~ltcabcuml)
lme.d<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl, data=jing.d, random=~ltcabcuml)
lme.e<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl, data=jing.e, random=~ltcabcuml)
lme.f<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl, data=jing.f, random=~ltcabcuml)



summary(lme.a)
summary(lme.b)
summary(lme.c)
summary(lme.d)
summary(lme.e)

##autocorrelation 

auto.a<-update(lme.a, correlation=corAR1())
auto.b<-update(lme.b, correlation=corAR1())
auto.c<-update(lme.c, correlation=corAR1())
auto.d<-update(lme.d, correlation=corAR1())
auto.e<-update(lme.e, correlation=corAR1())

summary(auto.a)
anova(auto.a, lme.a)
coef(auto.a, level=1)
plot(auto.a, idn~resid(.), abline=0)
plot(auto.a, resid(., type="p")~fitted(.)|idn, id=0.05, adj=-0.3)

summary(auto.a)
summary(auto.b)
summary(auto.c)
summary(auto.d)
summary(auto.e)

coef(auto.a)
coef(lme.a)


##interaction 
f<-read.dta("final1_before.dta")
g<-read.dta("final1_after.dta")
h<-read.dta("final2_before.dta")
i<-read.dta("final2_after.dta")
j<-read.dta("final3_before.dta")
k<-read.dta("final3_after.dta")
l<-read.dta("final4_before.dta")
m<-read.dta("final4_after.dta")
n<-read.dta("final5_before.dta")
o<-read.dta("final5_after.dta")

na.f <- !is.na( apply (f [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
f.nomiss <- f [na.f ,]
na.g <- !is.na( apply (g [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
g.nomiss <- g [na.g ,]
na.h <- !is.na( apply (h [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
h.nomiss <- h [na.h ,]
na.i <- !is.na( apply (i [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
i.nomiss <- i [na.i ,]
na.j <- !is.na( apply (j [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
j.nomiss <- j [na.j ,]
na.k <- !is.na( apply (k [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
k.nomiss <- k [na.k ,]
na.l <- !is.na( apply (l [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
l.nomiss <- l [na.l ,]
na.m <- !is.na( apply (m [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
m.nomiss <- m [na.m ,]
na.n <- !is.na( apply (n [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
n.nomiss <- n [na.n ,]
na.o <- !is.na( apply (o [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
o.nomiss <- o [na.o ,]


jing.f<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=f.nomiss)
jing.g<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=g.nomiss)
jing.h<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=h.nomiss)
jing.i<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=i.nomiss)
jing.j<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=j.nomiss)
jing.k<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=k.nomiss)
jing.l<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=l.nomiss)
jing.m<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=m.nomiss)
jing.n<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=n.nomiss)
jing.o<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=o.nomiss)


bef.f<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+ltcabcuml:social+ltcabcuml:conservative, data=jing.f, random=~ltcabcuml)
auto.f<-update(bef.f, correlation=corAR1())
summary(auto.f)
bef.g<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+ltcabcuml:social+ltcabcuml:conservative, data=jing.g, random=~ltcabcuml)
auto.g<-update(bef.g, correlation=corAR1())
summary(auto.g)
bef.h<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+ltcabcuml:social+ltcabcuml:conservative, data=jing.h, random=~ltcabcuml)
bef.i<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+ltcabcuml:social+ltcabcuml:conservative, data=jing.i, random=~ltcabcuml)
bef.j<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+ltcabcuml:social+ltcabcuml:conservative, data=jing.j, random=~ltcabcuml)
bef.k<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+ltcabcuml:social+ltcabcuml:conservative, data=jing.k, random=~ltcabcuml)
bef.l<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+ltcabcuml:social+ltcabcuml:conservative, data=jing.l, random=~ltcabcuml)
bef.m<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+ltcabcuml:social+ltcabcuml:conservative, data=jing.m, random=~ltcabcuml)
bef.n<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+ltcabcuml:social+ltcabcuml:conservative, data=jing.n, random=~ltcabcuml)
bef.o<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+ltcabcuml:social+ltcabcuml:conservative, data=jing.o, random=~ltcabcuml)

chr.g<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+cncabcuml:social+cncabcuml:conservative, data=jing.g, random=~cncabcuml)
summary(chr.g)

##before
summary(bef.f)
summary(bef.h)
summary(bef.j)
summary(bef.l)
summary(bef.n)

##after
summary(bef.g)
summary(bef.i)
summary(bef.k)
summary(bef.m)
summary(bef.o)



intauto.a<-update(lnt.a, correlation=corAR1())
summary(intauto.a)
anova(lnt.a, lme.a, test=F)
plot(lnt.a, idn~resid(.), abline=0)
plot(lme.a, idn~resid(.), abline=0)
plot(auto.a, idn~resid(.), abline=0)
plot(lnt.a, resid(., type="p")~fitted(.)|idn, id=0.05, adj=-0.3)

##no outliers
f<-read.dta("final1_noout.dta")
na.f <- !is.na( apply (f [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
f.nomiss <- f [na.f ,]
jing.f<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=f.nomiss)
lnt.f<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+ltcabcuml:social+ltcabcuml:conservative, data=jing.f, random=~ltcabcuml)
summary(lnt.f)

g<-read.dta("final1_before.dta")
na.g <- !is.na( apply (g [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
g.nomiss <- g [na.g ,]
jing.g<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=g.nomiss)
lnt.g<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+ltcabcuml:social+ltcabcuml:conservative, data=jing.g, random=~ltcabcuml)
summary(lnt.g)

g<-read.dta("final1_before.dta")
f<-read.dta("final1_after.dta")
na.f <- !is.na( apply (f [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
f.nomiss <- f [na.f ,]
jing.f<-groupedData(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl|idn, data=f.nomiss)
lnt.f<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl+social+conservative+ltcabcuml:social+ltcabcuml:conservative, data=jing.f, random=~ltcabcuml)
summary(lnt.f)

**autocorrelation 
ACF(auto.lme)
plot(ACF(auto.lme, maxLag=10), alpha=0.01)

auto1.lme<-update(auto.lme, correlation=corAR1())
auto1.lme
summary(auto1.lme)
anova(auto.lme, auto1.lme)

auto2.lme<-update(auto.lme, correlation=corARMA(q=2))
auto2.lme
anova(auto1.lme, auto2.lme, test=F)**AR1 better

auto3.lme<-update(auto.lme, correlation=corCAR1(form=~year))
anova(auto1.lme, auto3.lme, test=F)**same

auto4.lme<-update(auto.lme, correlation=corARMA(p=1, q=1))
auto4.lme
anova(auto1.lme, auto4.lme)**AR1 better
plot(ACF(auto4.lme, maxlag=10, resType="n"), alpha=0.01)**no sig AR after 1 lag)

**heteroskedasticity

hetero.lme<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl, data=jing.a, random=~ltcabcuml)
plot(hetero.lme, resid(.) ~ltcabcuml, abline=0)

hetero.lme<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl, data=jing.a, random=~ltcabcuml)
plot(hetero.lme, resid(.) ~ltcabcuml, abline=0)

decom.both<-lme(decom~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl, data=jing.a, random=~ltcabcuml+veto1l)
decom.both

uedecom<-lme(uescorecov~ltcabcuml+cncabcuml+veto1l+wmobcal + vturnl+ oldl+strikel+authlegl+lgdppcl+cpil+unempl+mill+dinvocl+openkl, data=jing.a, random=~ltcabcuml)
uedecom
summary(uedecom)
auto1.ue<-update(uedecom, correlation=corAR1())
auto1.ue
summary(auto1.ue)

######Pooling and No-Pooling Comparison 

##Complete pooling regression
lm.pooled1 <- lm (decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl, data=cws)
display(lm.pooled1)

##No pooling regression
country <- cws[,"idn"]

lm.unpooled1 <- lm (decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl + factor(country)-1, data=cws)
display (lm.unpooled1)

## Comparison between complete pooling and no-pooling 

y1 <- cws[,"decom"]
x1 <- cws[,"ltcabcuml"]
n1 <- length(y1)
x.jitter1 <- x1 + runif(n1,-.5,.5)

plot(x.jitter1,y1)
curve(coef(lm.pooled2)[1] + coef(lm.pooled2)[2]*x, add=TRUE)

display1 <- c(seq(1,8))

par(mfrow=c(2,4),mar=c(4,4,3,1),oma=c(1,1,2,1))
for (i in display1){
plot(x.jitter1[country==i],y1[country==i],xlim=c(0,50),ylim=c(0,50), xlab="left-cabinet",ylab="decom")
curve(coef(lm.pooled1)[1]+coef(lm.pooled1)[2]*x,lwd=.5,lty=2,add=TRUE)
curve(coef(lm.unpooled1)[i+14]+coef(lm.unpooled)[1]*x,lwd=.5,add=TRUE)
}

#######Varing Intercepts and Slopes (ALL)

y1a <- cws[,"decom"]
x1a <- cws[,"ltcabcuml"]
n1a <- length(y1)
country <- cws[,"idn"]
J1a <- 18

M31a <- lmer (decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl + (1 + x1 | country), data = cws)

display(M31a)

fixef(M31a)
ranef(M31a)


a.hat.M31a <- fixef(M31a)[1] + ranef(M31a)$country[,1]
b.hat.M31a <- fixef(M31a)[2] + ranef(M31a)$country[,2]

b.hat.unpooled.varying1a <- array (NA, c(J1a,15))
for (i in 1:J1a){
	lm.unpooled.varying1a <- lm(decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl,data=cws, subset=(country==i))
	b.hat.unpooled.varying1a[i,] <- coef(lm.unpooled.varying1a)
	}

lm.pooled1a <- lm (decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl,data=cws)

x.jitter1a <- x1a + runif(n1,-1,1)

##Social-democracy (ALL)

par (mfrow=c(2,3), mar=c(4,4,3,1), oma=c(1,1,2,1))
social.demo <- c (5,6,12,14,15)
for (j in social.demo){
  plot (x.jitter1a[country==j], y1a[country==j], xlim=c(-40,50), ylim=c(-40,50), xlab="left", ylab="decom")
  curve (coef(lm.pooled1a)[1] + coef(lm.pooled1a)[2]*x, lwd=.5, lty=2, col="gray10", add=TRUE)
  curve (b.hat.unpooled.varying1a[j,1] + b.hat.unpooled.varying1a[j,2]*x, lwd=.5, col="gray10", add=TRUE)
  curve (a.hat.M31a[j] + b.hat.M31a[j]*x,add=TRUE)
}


##Liberal (ALL)

par (mfrow=c(2,4), mar=c(4,4,3,1), oma=c(1,1,2,1))
liberal <- c (1, 4, 9, 11, 13, 16, 17, 18)
for (j in liberal){
  plot (x.jitter1a[country==j], y1a[country==j], xlim=c(-40,50), ylim=c(-40,50), xlab="left", ylab="decom")
  curve (coef(lm.pooled1a)[1] + coef(lm.pooled1a)[2]*x, lwd=.5, lty=2, col="gray10", add=TRUE)
  curve (b.hat.unpooled.varying1a[j,1] + b.hat.unpooled.varying1a[j,2]*x, lwd=.5, col="gray10", add=TRUE)
  curve (a.hat.M31a[j] + b.hat.M31a[j]*x,add=TRUE)
}


##Conservative (ALL)

par (mfrow=c(2,3), mar=c(4,4,3,1), oma=c(1,1,2,1))
conservative <- c (2, 3, 7, 8, 10)
for (j in conservative){
  plot (x.jitter1a[country==j], y1a[country==j], xlim=c(-40,50), ylim=c(-40,50), xlab="left", ylab="decom")
  curve (coef(lm.pooled1a)[1] + coef(lm.pooled1a)[2]*x, lwd=.5, lty=2, col="gray10", add=TRUE)
  curve (b.hat.unpooled.varying1a[j,1] + b.hat.unpooled.varying1a[j,2]*x, lwd=.5, col="gray10", add=TRUE)
  curve (a.hat.M31a[j] + b.hat.M31a[j]*x,add=TRUE)
}


##Estimation by regimes (ALL)
cws.regime <- read.dta("final1_regime.dta")

y11a <- cws.regime[,"decom"]
x11a <- cws.regime[,"ltcabcuml"]
n1a <- length(y11)
regime <- cws.regime[,"regime"]
J11a <- 3

x.jitter11a <- x11a + runif(n1,-1,1)

M311a <- lmer (decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl + (1 + x1 | regime), data = cws.regime)

display(M311a)

fixef(M311a)
ranef(M311a)

a.hat.M311a <- fixef(M311a)[1] + ranef(M311a)$regime[,1]
b.hat.M311a <- fixef(M311a)[2] + ranef(M311a)$regime[,2]

b.hat.unpooled.varying11a <- array (NA, c(J11a,15))
for (i in 1:J11a){
	lm.unpooled.varying11a <- lm(decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl,data=cws.regime, subset=(regime==i))
	b.hat.unpooled.varying11a[i,] <- coef(lm.unpooled.varying11a)
	}

lm.pooled11a <- lm (decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl,data=cws.regime)

x.jitter11a <- x11a + runif(n1,-1,1)


par (mfrow=c(1,3), mar=c(4,4,3,1), oma=c(1,1,2,1))
regime <- c (1,2,3)
for (j in regime){
  plot (x.jitter11a[regime==j], y11a[regime==j], xlim=c(-40,50), ylim=c(-40,50), xlab="left", ylab="decom")
  curve (coef(lm.pooled11a)[1] + coef(lm.pooled11a)[2]*x, lwd=.5, lty=2, col="gray10", add=TRUE)
  curve (b.hat.unpooled.varying11a[j,1] + b.hat.unpooled.varying11a[j,2]*x, lwd=.5, col="gray10", add=TRUE)
  curve (a.hat.M311a[j] + b.hat.M311a[j]*x, lwd=1,add=TRUE)
}


######Non-nested model by two criteria (Time & Regime)
##Regression for Time and Regime
##Control (All)
##plot$time(<1985)
time.temp <- read.dta("final1_before.dta")
na.ind.before <- !is.na( apply (time.temp [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
time <- time.temp [na.ind.before,]

yt1 <- time[,"decom"]
xt1 <- time[,"ltcabcuml"]
nt1 <- length(yt1)
country <- time[,"idn"]


MT1 <- lmer (decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl + social + conservative + ltcabcuml:social  + ltcabcuml:conservative + (1+ltcabcuml|idn), data=time)

display(MT1)
coef(MT1)
fixef(MT1)
ranef(MT1)

xt1.jittered <- xt1 + runif(nt1,-1,1)
a.hat.MT1 <- fixef(MT1)[1] + ranef(MT1)$idn[,1]
b.hat.MT1 <- fixef(MT1)[2] + ranef(MT1)$idn[,2]

plot(xt1.jittered,yt1,xlim=c(0,50),ylim=c(-50,50),xlab="Left Cabinet",ylab="Decommodation",main=paste("Overall:1985(-)"))
for (i in 1:18){
	curve(a.hat.MT1[i] + b.hat.MT1[i]*x, add=TRUE)
}

par (mfrow=c(2,3))

##Social Democracy

xt1.social <- c(xt1.jittered[country==5],xt1.jittered[country==6],xt1.jittered[country==12],xt1.jittered[country==14],xt1.jittered[country==15])
yt1.social <-
c(yt1[country==5],yt1[country==6],yt1[country==12],yt1[country==14],yt1[country==15])

plot(xt1.social,yt1.social,xlim=c(0,50),ylim=c(-50,50),xlab="Left Cabinet",ylab="Decommodation", main=paste("Social-Demo:1985(-)"))
curve(a.hat.MT1[5] + b.hat.MT1[5]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[6] + b.hat.MT1[6]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[12] + b.hat.MT1[12]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[14] + b.hat.MT1[14]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[15] + b.hat.MT1[15]*x, add=TRUE, lwd=1)

##Liberal

xt1.liberal <- c(xt1.jittered[country==1],xt1.jittered[country==4],xt1.jittered[country==9],xt1.jittered[country==11],xt1.jittered[country==13],xt1.jittered[country==16],xt1.jittered[country==17],xt1.jittered[country==18])
yt1.liberal <-
c(yt1[country==1],yt1[country==4],yt1[country==9],yt1[country==11],yt1[country==13],yt1[country==16],yt1[country==17],yt1[country==18])

plot(xt1.liberal,yt1.liberal,xlim=c(0,50),ylim=c(-50,50),xlab="Left Cabinet",ylab="Decommodation", main=paste("Liberal:1985(-)"))
curve(a.hat.MT1[1] + b.hat.MT1[1]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[4] + b.hat.MT1[4]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[9] + b.hat.MT1[9]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[11] + b.hat.MT1[11]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[13] + b.hat.MT1[13]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[16] + b.hat.MT1[16]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[17] + b.hat.MT1[17]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[18] + b.hat.MT1[18]*x, add=TRUE, lwd=1)

##Conservative
xt1.conservative <- c(xt1.jittered[country==2],xt1.jittered[country==3],xt1.jittered[country==7],xt1.jittered[country==8],xt1.jittered[country==10])
yt1.conservative <-
c(yt1[country==2],yt1[country==3],yt1[country==7],yt1[country==8],yt1[country==10])

plot(xt1.conservative,yt1.conservative,xlim=c(0,50),ylim=c(-50,50),xlab="Left Cabinet",ylab="Decommodation",main=paste("Conservative:1985(-)"))
curve(a.hat.MT1[2] + b.hat.MT1[2]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[3] + b.hat.MT1[3]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[7] + b.hat.MT1[7]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[8] + b.hat.MT1[8]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[10] + b.hat.MT1[10]*x, add=TRUE, lwd=1)

##plot$time(>1984)
time.temp <- read.dta("final1_after.dta")
na.ind.before <- !is.na( apply (time.temp [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
time <- time.temp [na.ind.before,]

yt1 <- time[,"decom"]
xt1 <- time[,"ltcabcuml"]
nt1 <- length(yt1)
country <- time[,"idn"]


MT1 <- lmer (decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl + social + conservative + ltcabcuml:social  + ltcabcuml:conservative + (1+ltcabcuml|idn), data=time)

display(MT1)
coef(MT1)
fixef(MT1)
ranef(MT1)

xt1.jittered <- xt1 + runif(nt1,-1,1)
a.hat.MT1 <- fixef(MT1)[1] + ranef(MT1)$idn[,1]
b.hat.MT1 <- fixef(MT1)[2] + ranef(MT1)$idn[,2]

##Social Democracy

xt1.social <- c(xt1.jittered[country==5],xt1.jittered[country==6],xt1.jittered[country==12],xt1.jittered[country==14],xt1.jittered[country==15])
yt1.social <-
c(yt1[country==5],yt1[country==6],yt1[country==12],yt1[country==14],yt1[country==15])

plot(xt1.social,yt1.social,xlim=c(0,50),ylim=c(-50,50),xlab="Left Cabinet",ylab="Decommodation", main=paste("Social-Demo:1985(+)"))
curve(a.hat.MT1[5] + b.hat.MT1[5]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[6] + b.hat.MT1[6]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[12] + b.hat.MT1[12]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[14] + b.hat.MT1[14]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[15] + b.hat.MT1[15]*x, add=TRUE, lwd=1)

##Liberal

xt1.liberal <- c(xt1.jittered[country==1],xt1.jittered[country==4],xt1.jittered[country==9],xt1.jittered[country==11],xt1.jittered[country==13],xt1.jittered[country==16],xt1.jittered[country==17],xt1.jittered[country==18])
yt1.liberal <-
c(yt1[country==1],yt1[country==4],yt1[country==9],yt1[country==11],yt1[country==13],yt1[country==16],yt1[country==17],yt1[country==18])

plot(xt1.liberal,yt1.liberal,xlim=c(0,50),ylim=c(-50,50),xlab="Left Cabinet",ylab="Decommodation", main=paste("Liberal:1985(+)"))
curve(a.hat.MT1[1] + b.hat.MT1[1]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[4] + b.hat.MT1[4]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[9] + b.hat.MT1[9]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[11] + b.hat.MT1[11]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[13] + b.hat.MT1[13]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[16] + b.hat.MT1[16]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[17] + b.hat.MT1[17]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[18] + b.hat.MT1[18]*x, add=TRUE, lwd=1)

##Conservative
xt1.conservative <- c(xt1.jittered[country==2],xt1.jittered[country==3],xt1.jittered[country==7],xt1.jittered[country==8],xt1.jittered[country==10])
yt1.conservative <-
c(yt1[country==2],yt1[country==3],yt1[country==7],yt1[country==8],yt1[country==10])

plot(xt1.conservative,yt1.conservative,xlim=c(0,50),ylim=c(-50,50),xlab="Left Cabinet",ylab="Decommodation",main=paste("Conservative:1985(+)"))
curve(a.hat.MT1[2] + b.hat.MT1[2]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[3] + b.hat.MT1[3]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[7] + b.hat.MT1[7]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[8] + b.hat.MT1[8]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[10] + b.hat.MT1[10]*x, add=TRUE, lwd=1)

#### Overall:1985(+)
plot(xt1.jittered,yt1,xlim=c(0,50),ylim=c(-50,50),xlab="Left Cabinet",ylab="Decommodation",main=paste("Overall:1985(+)"))
for (i in 1:18){
	curve(a.hat.MT1[i] + b.hat.MT1[i]*x, add=TRUE)
}

######Plot$before&after by regimes

par(mfrow=c(1,2))
##before 1985
time.temp <- read.dta("final1_before.dta")
na.ind.before <- !is.na( apply (time.temp [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
time <- time.temp [na.ind.before,]

yt1 <- time[,"decom"]
xt1 <- time[,"ltcabcuml"]
nt1 <- length(yt1)
country <- time[,"idn"]

MT1 <- lmer (decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl + (1+ltcabcuml|regime), data=time)

display(MT1)
coef(MT1)
fixef(MT1)
ranef(MT1)

xt1.jittered <- xt1 + runif(nt1,-1,1)
a.hat.MT1 <- fixef(MT1)[1] + ranef(MT1)$regime[,1]
b.hat.MT1 <- fixef(MT1)[2] + ranef(MT1)$regime[,2]

plot(xt1.jittered,yt1,xlim=c(0,50),ylim=c(-20,100),xlab="left-cabinet",ylab="Decommodification")
curve(a.hat.MT1[1] + b.hat.MT1[1]*x, add=TRUE)
curve(a.hat.MT1[2] + b.hat.MT1[2]*x, add=TRUE,lty=3)
curve(a.hat.MT1[3] + b.hat.MT1[3]*x, add=TRUE,lty=4)

##After1985

time.temp <- read.dta("final1_after.dta")
na.ind.before <- !is.na( apply (time.temp [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
time <- time.temp [na.ind.before,]

yt1 <- time[,"decom"]
xt1 <- time[,"ltcabcuml"]
nt1 <- length(yt1)
country <- time[,"idn"]

MT1 <- lmer (decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl + (1+ltcabcuml|regime), data=time)

display(MT1)
coef(MT1)
fixef(MT1)
ranef(MT1)

xt1.jittered <- xt1 + runif(nt1,-1,1)
a.hat.MT1 <- fixef(MT1)[1] + ranef(MT1)$regime[,1]
b.hat.MT1 <- fixef(MT1)[2] + ranef(MT1)$regime[,2]

plot(xt1.jittered,yt1,xlim=c(0,50),ylim=c(-20,100),xlab="left-cabinet",ylab="Decommodification")
curve(a.hat.MT1[1] + b.hat.MT1[1]*x, add=TRUE)
curve(a.hat.MT1[2] + b.hat.MT1[2]*x, add=TRUE,lty=3)
curve(a.hat.MT1[3] + b.hat.MT1[3]*x, add=TRUE,lty=4)


##No-Control

par(mfrow=c(1,2))

time.temp <- read.dta("final1_before.dta")
na.ind.before <- !is.na( apply (time.temp [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
time <- time.temp [na.ind.before,]

yt1 <- time[,"decom"]
xt1 <- time[,"ltcabcuml"]
nt1 <- length(yt1)
country <- time[,"idn"]

MT1 <- lmer (decom ~ ltcabcuml + (1+ltcabcuml|regime), data=time)

display(MT1)
coef(MT1)
fixef(MT1)
ranef(MT1)

xt1.jittered <- xt1 + runif(nt1,-1,1)
a.hat.MT1 <- fixef(MT1)[1] + ranef(MT1)$regime[,1]
b.hat.MT1 <- fixef(MT1)[2] + ranef(MT1)$regime[,2]

plot(xt1.jittered,yt1,xlim=c(0,50),ylim=c(-20,100),xlab="left-cabinet",ylab="Decommodification")
curve(a.hat.MT1[1] + b.hat.MT1[1]*x, add=TRUE)
curve(a.hat.MT1[2] + b.hat.MT1[2]*x, add=TRUE,lty=3)
curve(a.hat.MT1[3] + b.hat.MT1[3]*x, add=TRUE,lty=4)


time.temp <- read.dta("final1_after.dta")
na.ind.before <- !is.na( apply (time.temp [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
time <- time.temp [na.ind.before,]

yt1 <- time[,"decom"]
xt1 <- time[,"ltcabcuml"]
nt1 <- length(yt1)
country <- time[,"idn"]

MT1 <- lmer (decom ~ ltcabcuml + (1+ltcabcuml|regime), data=time)

display(MT1)
coef(MT1)
fixef(MT1)
ranef(MT1)

xt1.jittered <- xt1 + runif(nt1,-1,1)
a.hat.MT1 <- fixef(MT1)[1] + ranef(MT1)$regime[,1]
b.hat.MT1 <- fixef(MT1)[2] + ranef(MT1)$regime[,2]

plot(xt1.jittered,yt1,xlim=c(0,50),ylim=c(-20,100),xlab="left-cabinet",ylab="Decommodification")
curve(a.hat.MT1[1] + b.hat.MT1[1]*x, add=TRUE)
curve(a.hat.MT1[2] + b.hat.MT1[2]*x, add=TRUE,lty=3)
curve(a.hat.MT1[3] + b.hat.MT1[3]*x, add=TRUE,lty=4)

######Non-nested model by controlling time and countries 

time.temp <- read.dta("final1_regime.dta")
na.ind.before <- !is.na( apply (time.temp [,c("decom", "ltcabcuml", "cncabcuml", "veto1l","wmobcal", "vturnl", "oldl", "strikel", "authlegl", "lgdppcl", "cpil", "unempl", "mill", "dinvocl", "openkl")],1, sum))
time <- time.temp [na.ind.before,]

MT1 <- lmer (decom ~ ltcabcuml + cncabcuml + veto1l + wmobcal + vturnl + oldl + strikel + authlegl + lgdppcl + cpil + unempl + mill + dinvocl + openkl + (1+ltcabcuml|idn), data=time1)

head(time)
time[1:2] <- lapply(time[1:2],factor)

display(MT1)
coef(MT1)
fixef(MT1)
ranef(MT1)

xt1.jittered <- xt1 + runif(nt1,-1,1)
a.hat.MT1 <- fixef(MT1)[1] + ranef(MT1)$idn[,1]
b.hat.MT1 <- fixef(MT1)[2] + ranef(MT1)$idn[,2]

plot(xt1.jittered,yt1,xlim=c(0,50),ylim=c(-50,50),xlab="Left Cabinet",ylab="Decommodation",main=paste("Overall:1985(-)"))
for (i in 1:18){
	curve(a.hat.MT1[i] + b.hat.MT1[i]*x, add=TRUE)
}

par (mfrow=c(2,3))

##Social Democracy

xt1.social <- c(xt1.jittered[country==5],xt1.jittered[country==6],xt1.jittered[country==12],xt1.jittered[country==14],xt1.jittered[country==15])
yt1.social <-
c(yt1[country==5],yt1[country==6],yt1[country==12],yt1[country==14],yt1[country==15])

plot(xt1.social,yt1.social,xlim=c(0,50),ylim=c(-50,50),xlab="Left Cabinet",ylab="Decommodation", main=paste("Social-Demo:1985(-)"))
curve(a.hat.MT1[5] + b.hat.MT1[5]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[6] + b.hat.MT1[6]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[12] + b.hat.MT1[12]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[14] + b.hat.MT1[14]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[15] + b.hat.MT1[15]*x, add=TRUE, lwd=1)

##Liberal

xt1.liberal <- c(xt1.jittered[country==1],xt1.jittered[country==4],xt1.jittered[country==9],xt1.jittered[country==11],xt1.jittered[country==13],xt1.jittered[country==16],xt1.jittered[country==17],xt1.jittered[country==18])
yt1.liberal <-
c(yt1[country==1],yt1[country==4],yt1[country==9],yt1[country==11],yt1[country==13],yt1[country==16],yt1[country==17],yt1[country==18])

plot(xt1.liberal,yt1.liberal,xlim=c(0,50),ylim=c(-50,50),xlab="Left Cabinet",ylab="Decommodation", main=paste("Liberal:1985(-)"))
curve(a.hat.MT1[1] + b.hat.MT1[1]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[4] + b.hat.MT1[4]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[9] + b.hat.MT1[9]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[11] + b.hat.MT1[11]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[13] + b.hat.MT1[13]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[16] + b.hat.MT1[16]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[17] + b.hat.MT1[17]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[18] + b.hat.MT1[18]*x, add=TRUE, lwd=1)

##Conservative
xt1.conservative <- c(xt1.jittered[country==2],xt1.jittered[country==3],xt1.jittered[country==7],xt1.jittered[country==8],xt1.jittered[country==10])
yt1.conservative <-
c(yt1[country==2],yt1[country==3],yt1[country==7],yt1[country==8],yt1[country==10])

plot(xt1.conservative,yt1.conservative,xlim=c(0,50),ylim=c(-50,50),xlab="Left Cabinet",ylab="Decommodation",main=paste("Conservative:1985(-)"))
curve(a.hat.MT1[2] + b.hat.MT1[2]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[3] + b.hat.MT1[3]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[7] + b.hat.MT1[7]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[8] + b.hat.MT1[8]*x, add=TRUE, lwd=1)
curve(a.hat.MT1[10] + b.hat.MT1[10]*x, add=TRUE, lwd=1)
